% I've (fede) modified this function and eliminated redundancies and
% adapted to Aiyagari files

function  [alpha_star, mu_t0, polhat, eta_t, sigma_etasq, mu0_ya_hist, hist_mat,Linv,mu_t0_hist] = bc_update_iid_fast(r, W_cc, kappa, y_t, a_t, sigma_c,psi, cov_fun, ex_signals, eps_eta,b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c, psi_tau)

%  inputs
%  parameters: r, ybar, rho_y, W_cc, kappa
%  y_t = current state y_t
%  a_t = current state a_t
%  y_hist = past states y_1 ~ y_(t-1)
%  a_hist = past states a_1 ~ a_(t-1)
%  eta_hist = past signals eta_1 ~ eta_(t-1)
%  sigma_etasq_hist = past signal precision^(-1), after forgetting
%  prior functions and parameters: mu0, sigma_c, psi, cov_fun, ex_signals, res_action
%  

% NOTE: mu0 must be the time-zero prior policy function as well as the
% "correct" policy function

%  1) Beliefs based on eta_{t-1}

%  Ignore any point with no information
%y_hist = y_hist(sigma_etasq_hist < 1e10);
%a_hist = a_hist(sigma_etasq_hist < 1e10);
%eta_hist = eta_hist(sigma_etasq_hist < 1e10);
%sigma_etasq_hist = sigma_etasq_hist(sigma_etasq_hist < 1e10);



mu0_ytat = fast_interp( y_t + (1+r)*a_t, ctilde, stateGrid_step, stateGrid,stateGrid_size);

%  prior on chat(y) and Sigma(y,y')
if isempty(hist_mat)
    
    Sigma_t0 = sigma_c^2;     %Var(cstar(y) | eta_hist )
    mu_t0    = mu0_ytat;
    
else
    
    y_hist   = hist_mat(:,1);
    a_hist   = hist_mat(:,2);
    eta_hist = hist_mat(:,3);
    sigma_etasq_hist = hist_mat(:,4);
    
    if cov_fun == 1
        
        temp1_yat =  (1+r)*a_t + y_t - (1+r)*a_hist' - y_hist';
        Sigma12   = sigma_c^2*exp(-psi*(temp1_yat.^2) - psi_tau);
        
        
        v = Linv*Sigma12';
        sum_vsq = sum(v.^2);
        Sigma_t0 = sigma_c^2 - sum_vsq;
        
        mu_t0 = mu0_ytat + Sigma12*mu_t0_hist;
        
    elseif cov_fun ==2 
        temp1_yat = bsxfun(@minus, (1+r)*a_t + y_t, (1+r)*a_hist' + y_hist');
        Sigma12  = sigma_c^2*exp(-psi*(abs(temp1_yat)));
        
        
        v = Linv*Sigma12';
        sum_vsq = sum(v.^2);
        Sigma_t0 = sigma_c^2 - sum_vsq;   %Var(cstar(y) | eta_hist )
        
        mu_t0 = mu0_ytat + Sigma12*mu_t0_hist;

    end
    
end



%  2) updating with current signal eta_t

%  applying the optimal signal variance \sigma_{\eta,t}^2 the resulting signal-to-noise ratio is:

%  alpha_star = max((Sigma_t0 - kappa/(W_cc*2))./Sigma_t0,0);

if ex_signals > 0
    sigma_etasq = ex_signals;
    alpha_star  = Sigma_t0/(Sigma_t0 + sigma_etasq);
    
    eta_t = mu0_ytat + sqrt(sigma_etasq)*eps_eta;

    mu0_ya_hist = [mu0_ya_hist; mu0_ytat];
    if isempty(hist_mat)
        y_hist = y_t;
        a_hist = a_t;
        eta_hist = eta_t;
        sigma_etasq_hist = sigma_etasq;
        
        Linv = 1/sqrt( sigma_c^2 + sigma_etasq);
        
    else
        
        y_hist = [y_hist; y_t];
        a_hist = [a_hist; a_t];
        eta_hist = [eta_hist;eta_t];
        sigma_etasq_hist = [sigma_etasq_hist; sigma_etasq];
        
        
        D    = sqrt( sigma_c^2 + sigma_etasq - sum_vsq);
        temp = Linv;
        n    = length(y_hist);
        
        Linv(n, n) = 1/D;
        Linv(n, 1:n-1) = -v'*temp/D;
    end
    
    hist_mat = [ y_hist, a_hist, eta_hist, sigma_etasq_hist];
    mu_t0_hist = (Linv')*(Linv*(eta_hist - mu0_ya_hist));
    
elseif Sigma_t0 > kappa/(2*W_cc)
    sigma_etasq = (kappa/(2*W_cc))*Sigma_t0/(Sigma_t0 - kappa/(2*W_cc));
    alpha_star = max((Sigma_t0 - kappa/(W_cc*2))./Sigma_t0,0);
    
    eta_t = mu0_ytat + sqrt(sigma_etasq)*eps_eta;
    
    
    mu0_ya_hist = [mu0_ya_hist; mu0_ytat];
    if isempty(hist_mat)
        y_hist = y_t;
        a_hist = a_t;
        eta_hist = eta_t;
        sigma_etasq_hist = sigma_etasq;
        
        Linv = 1/sqrt( sigma_c^2 + sigma_etasq);
        
    else
        
        y_hist = [y_hist; y_t];
        a_hist = [a_hist; a_t];
        eta_hist = [eta_hist;eta_t];
        sigma_etasq_hist = [sigma_etasq_hist; sigma_etasq];
        
        
        D    = sqrt( sigma_c^2 + sigma_etasq - sum_vsq);
        temp = Linv;
        n    = length(y_hist);
        
        Linv(n, n) = 1/D;
        Linv(n, 1:n-1) = -v'*temp/D;
    end
    
    hist_mat = [ y_hist, a_hist, eta_hist, sigma_etasq_hist];
    mu_t0_hist = (Linv')*(Linv*(eta_hist - mu0_ya_hist));
    
else
    sigma_etasq = 1e10;
    alpha_star = max((Sigma_t0 - kappa/(W_cc*2))./Sigma_t0,0);
    eps_eta = 0;
    
    eta_t = mu0_ytat + sqrt(sigma_etasq)*eps_eta;
    
end


polhat = mu_t0.*(1-alpha_star) + alpha_star.*eta_t;

% polhat is the choice of c. if this violates the borrowing constraint then
% the agent consumes as much as possible
if log_c == 0
    if  polhat > (1+r).*a_t + y_t - b
        polhat = (1+r).*a_t + y_t - b;
    end
    
    if  polhat < 0
        polhat = 1e-2;
    end
    
    
    
else
    if polhat > log((1+r).*a_t + y_t - b)
        polhat = log((1+r).*a_t + y_t - b);
    end
    
    
end

end